home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / pascal / tpl60n19.zip / TESTPRGS.ZIP / DLOG.PAS < prev    next >
Pascal/Delphi Source File  |  1993-02-14  |  10KB  |  352 lines
  1. PROGRAM DLog;   { ported from Fortran original 05-01-92 Norbert Juffa }
  2.  
  3. {$A+,B-,D-,E+,F-,G-,I-,L-,N-,O-,R-,S-,V-,X-}
  4.  
  5. USES MachArit;
  6.  
  7. {
  8. C     PROGRAM TO TEST DLOG
  9. C
  10. C     DATA REQUIRED
  11. C
  12. C        NONE
  13. C
  14. C     SUBPROGRAMS REQUIRED FROM THIS PACKAGE
  15. C
  16. C        MACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING
  17. C                 INFORMATION ON THE FLOATING-POINT ARITHMETIC
  18. C                 SYSTEM.  NOTE THAT THE CALL TO MACHAR CAN
  19. C                 BE DELETED PROVIDED THE FOLLOWING FOUR
  20. C                 PARAMETERS ARE ASSIGNED THE VALUES INDICATED
  21. C
  22. C                 IBETA - THE RADIX OF THE FLOATING-POINT SYSTEM
  23. C                 IT    - THE NUMBER OF BASE-IBETA DIGITS IN THE
  24. C                         SIGNIFICAND OF A FLOATING-POINT NUMBER
  25. C                 XMIN  - THE SMALLEST NON-VANISHING FLOATING-POINT
  26. C                         POWER OF THE RADIX
  27. C                 XMAX  - THE LARGEST FINITE FLOATING-POINT NO.
  28. C
  29. C        REN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL
  30. C                 NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1)
  31. C
  32. C
  33. C     STANDARD FORTRAN SUBPROGRAMS REQUIRED
  34. C
  35. C         DABS, DLOG, DLOG10, DMAX1, DFLOAT, DSIGN, DSQRT
  36. C
  37. C
  38. C     LATEST REVISION - DECEMBER 6, 1979
  39. C
  40. C     AUTHOR - W. J. CODY
  41. C              ARGONNE NATIONAL LABORATORY
  42. C
  43. C
  44. }
  45.  
  46.  
  47.  
  48. FUNCTION REN (K: LONGINT): REAL;
  49.  
  50. {
  51.       DOUBLE PRECISION FUNCTION REN(K)
  52. C
  53. C     RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 266 BY PIKE AND
  54. C      HILL (MODIFIED BY HANSSON), COMMUNICATIONS OF THE ACM,
  55. C      VOL. 8, NO. 10, OCTOBER 1965.
  56. C
  57. C     THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTERS WITH
  58. C      FIXED POINT WORDLENGTH OF AT LEAST 29 BITS.  IT IS
  59. C      BEST IF THE FLOATING POINT SIGNIFICAND HAS AT MOST
  60. C      29 BITS.
  61. C
  62. }
  63.  
  64. VAR   J:  LONGINT;
  65. CONST IY: LONGINT = 100001;
  66.  
  67. BEGIN
  68.    J  := K;
  69.    IY := IY * 125;
  70.    IY := IY - (IY DIV 2796203) * 2796203;
  71.    REN:= 1.0 * (IY) / 2796203.0e0 * (1.0e0 + 1.0e-6 + 1.0e-12);
  72. END;
  73.  
  74.  
  75. FUNCTION LOG (X: REAL): REAL;
  76. BEGIN
  77.    LOG  := LN (X) * 0.43429448190325182765112891891660508;
  78. END;
  79.  
  80.  
  81. FUNCTION MAX1 (A, B:REAL): REAL;
  82. BEGIN
  83.    IF A > B THEN
  84.       MAX1 := A
  85.    ELSE
  86.       MAX1 := B;
  87. END;
  88.  
  89.  
  90.  
  91. VAR   I,IBETA,IEXP,IOUT,IRND,IT,I1,J,K1,K2,
  92.       K3,MACHEP,MAXEXP,MINEXP,N,NEGEP,NGRD: LONGINT;
  93.  
  94.       A,AIT,ALBETA,B,BETA,C,DEL,EIGHT,EPS,
  95.       EPSNEG,HALF,ONE,T, R6,R7,TENTH,W,X,
  96.       XL,XMAX,XMIN,XN,X1,Y,Z,ZERO,ZZ,FOUR,
  97.       TWO,THREE,NINETENTH,FIFTEEN,SIXTEEN,
  98.       TWENTYONE,THIRTYONE,TWOHUNDREDFORTY,
  99.       FIVEHUNDREDTWELVE:                    REAL;
  100.  
  101. LABEL 100, 110, 120, 150, 160, 220, 230, 240, 300;
  102.  
  103. BEGIN
  104.  
  105.    N := 1000000;  { number of trials }
  106.  
  107.    MACHAR (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP,
  108.            EPS,EPSNEG,XMIN,XMAX);
  109.    PRINTPARAM (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP,
  110.                EPS,EPSNEG,XMIN,XMAX);
  111.    BETA      := IBETA;
  112.    ALBETA    := LN (BETA);
  113.    AIT       := IT;
  114.    J         := IT DIV 3;
  115.    ZERO      := 0;
  116.    ONE       := 1;
  117.    TWO       := 2;
  118.    THREE     := 3;
  119.    FOUR      := 4;
  120.    EIGHT     := 8;
  121.    FIFTEEN   := 15;
  122.    SIXTEEN   := 16;
  123.    THIRTYONE := 31;
  124.    TWENTYONE := 21;
  125.    TENTH     := 0.1;
  126.    HALF      := 0.5;
  127.    NINETENTH := 0.9;
  128.    TWOHUNDREDFORTY  := 240;
  129.    FIVEHUNDREDTWELVE:= 512;
  130.    C         := ONE;
  131.  
  132.    FOR I := 1 TO J DO BEGIN
  133.       C := C / BETA;
  134.    END;
  135.  
  136.    B  := ONE + C;
  137.    A  := ONE - C;
  138.    XN := N;
  139.    I1 := 0;
  140.  
  141. {-----------------------------------------------------------------}
  142. {      RANDOM ARGUMENT ACCURACY TESTS                             }
  143. {-----------------------------------------------------------------}
  144.  
  145.    FOR J := 1 TO 4 DO BEGIN
  146.       K1 := 0;
  147.       K3 := 0;
  148.       X1 := ZERO;
  149.       R6 := ZERO;
  150.       R7 := ZERO;
  151.       DEL:= (B - A) / XN;
  152.       XL := A;
  153.  
  154.       FOR I := 1 TO N DO BEGIN
  155.          X  := DEL * REN(I1) + XL;
  156.          IF J <> 1 THEN
  157.             GOTO 100;
  158.          Y := X - HALF;
  159.          Y := Y - HALF;
  160.          ZZ:= LN (X);
  161.          Z := (Y * (ONE / THREE - Y / FOUR) - HALF) * Y * Y + Y;
  162.          GOTO 150;
  163. 100:     IF J <> 2 THEN
  164.             GOTO 110;
  165.          X := X + EIGHT;
  166.          X := X - EIGHT;
  167.          Y := X / SIXTEEN;
  168.          Y := X + Y;
  169.          Z := LN (X);
  170.          ZZ:= LN (Y);
  171.          ZZ:= ZZ - 7.7746816434842581e-5; { Ln (17/16) - 31/512) }
  172.          ZZ:= ZZ - THIRTYONE/FIVEHUNDREDTWELVE;
  173.          GOTO 150;
  174. 110:     IF J <> 3 THEN
  175.             GOTO 120;
  176.          X := X + EIGHT;
  177.          X := X - EIGHT;
  178.          T := X * TENTH;
  179.          Y := X + T;
  180.          Z := LOG (X);
  181.          ZZ:= LOG (Y);
  182.          ZZ:= ZZ - 3.7706015822504075e-4;  { Log10 (11/10) - 21/512) }
  183.          ZZ:= ZZ - TWENTYONE/FIVEHUNDREDTWELVE;
  184.          GOTO 150;
  185. 120:     T := X * X;
  186.          Z := LN (T);
  187.          ZZ:= LN (X);
  188.          ZZ:= ZZ + ZZ;
  189.  
  190. 150:     IF Z <> ZERO THEN
  191.             W := (Z - ZZ) / Z
  192.          ELSE IF ZZ <> ZERO THEN
  193.             W := ONE;
  194.          IF W > ZERO THEN
  195.             K1 := K1 + 1;
  196.          IF W < ZERO THEN
  197.             K3 := K3 + 1;
  198.          W := ABS (W);
  199.          IF W <= R6 THEN
  200.             GOTO 160;
  201.          R6 := W;
  202.          X1 := X;
  203. 160:     R7 := R7 + W * W;
  204.          XL := XL + DEL;
  205.       END;
  206.  
  207.       K2 := N - K3 - K1;
  208.       R7 := SQRT (R7/XN);
  209.  
  210.       IF J = 1 THEN BEGIN
  211.          WRITELN;
  212.          WRITELN ;
  213.          WRITELN ('TEST OF LN (X) VS T.S. EXPANSION OF LN(1+Y)');
  214.          WRITELN;
  215.       END;
  216.       IF J = 2 THEN BEGIN
  217.          WRITELN;
  218.          WRITELN;
  219.          WRITELN ('TEST OF LN(X) VS LN(17X/16)-LN(17/16)');
  220.          WRITELN;
  221.       END;
  222.       IF J = 3 THEN BEGIN
  223.          WRITELN;
  224.          WRITELN;
  225.          WRITELN ('TEST OF LOG10(X) VS LOG10(11X/10)-LOG10(11/10)');
  226.          WRITELN;
  227.          END;
  228.       IF J = 4 THEN BEGIN
  229.          WRITELN;
  230.          WRITELN;
  231.          WRITELN ('TEST OF LN (X*X) VS 2*LN(X)');
  232.          WRITELN;
  233.       END;
  234.       IF J = 1 THEN BEGIN
  235.          WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL');
  236.          WRITELN ('(1-EPS,1+EPS), WHERE EPS = ', C);
  237.          WRITELN;
  238.       END;
  239.       IF J <> 1 THEN BEGIN
  240.          WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL');
  241.          WRITELN ('(', A, ',', B, ')');
  242.          WRITELN;
  243.       END;
  244.       IF J <> 3 THEN BEGIN
  245.          WRITELN ('LN (X) WAS LARGER', K1:6, ' TIMES');
  246.          WRITELN ('           AGREED', K2:6, ' TIMES');
  247.          WRITELN ('  AND WAS SMALLER', K3:6, ' TIMES');
  248.       END;
  249.       IF J = 3 THEN BEGIN
  250.          WRITELN ('LOG (X) WAS LARGER', K1:6, ' TIMES');
  251.          WRITELN ('            AGREED', K2:6, ' TIMES');
  252.          WRITELN ('   AND WAS SMALLER', K3:6, ' TIMES');
  253.       END;
  254.       WRITELN;
  255.       WRITELN ('THERE ARE ', IT, ' BASE ', IBETA,
  256.                ' SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER');
  257.       WRITELN;
  258.       W := -999;
  259.       IF R6 <> ZERO THEN
  260.          W := LN (ABS(R6))/ALBETA;
  261.       WRITELN ('THE MAXIMUM RELATIVE ERROR OF          ', R6:12,
  262.                ' = ', IBETA, ' **', W:7:2);
  263.       WRITELN ('OCCURED FOR X = ', X1);
  264.       W := MAX1 (AIT+W,ZERO);
  265.       WRITELN;
  266.       WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
  267.                ' SIGNIFICANT DIGITS IS        ', W:7:2);
  268.       W := -999.0;
  269.       IF R7 <> ZERO THEN
  270.          W := LN (ABS(R7))/ALBETA;
  271.       WRITELN;
  272.       WRITELN ('THE ROOT MEAN SQUARE RELATIVE ERROR WAS', R7:12,
  273.                ' = ', IBETA, ' **' , W:7:2);
  274.       W := MAX1 (AIT+W,ZERO);
  275.       WRITELN;
  276.       WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
  277.                ' SIGNIFICANT DIGITS IS        ', W:7:2);
  278.       IF J > 1 THEN
  279.          GOTO 230;
  280.       A := SQRT (HALF);
  281.       B := FIFTEEN / SIXTEEN;
  282.       GOTO 300;
  283. 230:  IF J > 2 THEN
  284.          GOTO 240;
  285.       A := SQRT (TENTH);
  286.       B := NINETENTH;
  287.       GOTO 300;
  288. 240:  A := SIXTEEN;
  289.       B := TWOHUNDREDFORTY;
  290. 300:
  291.    END;
  292.  
  293. {-----------------------------------------------------------------}
  294. {      SPECIAL TESTS                                              }
  295. {-----------------------------------------------------------------}
  296.  
  297.    WRITEL